arXiv:hep-ph/0512210v2 17 Aug 2006 


hep-ph/0512210 
LPTHE-05-32 
December 2005 
revised August 2006 


Dispelling the myth for the kt jet-finder 

Matteo Cacciari and Gavin P. Salam 

LPTHE, Universities of Paris VI and VII and CNRS, 

Abstract 

At high-energy colliders, jets of hadrons are the observable counterparts of the 
perturbative concepts of quarks and gluons. Good procedures for identifying jets 
are central to experimental analyses and comparisons with theory. The kt family of 
successive recombination jet finders has been widely advocated because of its con¬ 
ceptual simplicity and flexibility and its unique ability to approximately reconstruct 
the partonic branching sequence in an event. Until now however, it had been be¬ 
lieved that for an ensemble of N particles the algorithmic complexity of the kt jet 
finder scaled as a severe issue in the high multiplicity environments of LHC 
and heavy-ion colliders. We here show that the computationally complex part of 
kt jet-clustering can be reduced to two-dimensional nearest neighbour location for a 
dynamic set of points. Borrowing techniques developed for this extensively studied 
problem in computational geometry, kt jet-finding can then be performed in Ain A 
time. Code based on these ideas is found to run faster than all other jet finders in 
current use. 


1 Introduction 

Partons (quarks and gluons), are the concepts that are central to discussions of the QCD 
aspects of high-energy collisions such as those at the Fermilab Tevatron and the future 
Large Hadron Collider (LHC) at CERN. Quarks and gluons, however, are not observable, 
and in their place one sees jets, collimated bunches of high-energy hadrons which are the 
result of the fragmentation and hadronisation of the original hard (high-energy) partons. 
Today’s limited understanding of non-perturbative QCD is such that it is not currently 
possible to predict the exact patterns of hadrons produced. Instead one makes predictions 
in terms of quarks and gluons and relates these to observations in terms of hadron jets. 



Naively, jets are easily identified — one simply searches for bunches of collimated 
hadrons. However, to carry out accurate comparisons between parton-level predictions and 
hadron-level observations one needs a well-dehned ‘jet-hnding’ procedure. The jet-hnder is 
applied both to perturbatively predicted partonic conhgurations and to observed hadronic 
conhgurations and one then directly compares distributions for the predicted partonic jets 
and the observed hadronic jets. Though partonic and hadronic jets are not equivalent, 
there is strong evidence (theoretical [1] and experimental [2]) that the comparison can be 
performed with controlled accuracy. 

Insofar as jet-hnding is an approximate attempt to invert the quantum mechanical 
processes of QCD branching and hadronisation, it is not a unique procedure. Various kinds 
of jet-hnders have been proposed, among them cone-type [1,3] and sequential-clustering 
[4-7] jet-hnders (for alternatives, see [8-11]). 

Cone jet-hnders are the most frequently used at the Tevatron. They are based on 
identifying energy-how into cones in (pseudo)rapidity rj = — lntan0/2 and azimuth 0, 
together with various steps of iteration, merging and splitting of the cones to obtain the 
hnal jets. Cone jet-hnders tend to be rather complex, diherent experiments have used 
diherent variants (some of them infrared unsafe), and it is often difficult to know exactly 
which jet-hnder to use in theoretical comparisons. 

In contrast, the cluster-type jet-hnders, generally based on successive pair-wise recombi¬ 
nation of particles, have simple dehnitions and are all infrared safe (for reviews see [12,13]). 
We shall focus here on the most widely used of them, the kt jet-hnder [5], dehned below. 
Among its physics advantages are (a) that it purposely mimics a walk backwards through 
the QCD branching sequence, which means that reconstructed jets naturally collect most 
of the particles radiated from an original hard parton, giving better particle mass mea¬ 
surements [14,15], general kinematic reconstruction [16] and gaps-between-jets identihca- 
tion [17] (of relevance to Higgs searches); and (b) it allows one to decompose a jet into 
constituent subjets, which is useful for identifying decay products of fast-moving heavy 
particles (see e.g. [18]) and various QCD studies. This has led to the widespread adoption 
of the kt jet-hnder in the LEP (e’''e“ collisions) and HERA (ep) communities. 

Despite its advantages, kt clustering has so far seen only limited study [19-21] at the 
Tevatron. The reasons for this are not entirely clear. One known drawback of the kt 
jet hnder for high-multiplicity hadron-collider environments is its apparent algorithmic 
slowness: to cluster N particles into jets requires O {N^) operations in current implemen¬ 
tations [22]. For a typical event at the upcoming LHC, with an expected multiplicity of 
N = O (2000), this translates into a clustering time of O (10 s) of CPU time on a modern 
O (3 GHz) processor; this is considerable given that the clustering has to be repeated for 
millions of events. For a typical heavy-ion event at LHC, where N = O (50000), the clus¬ 
tering time would grow to an unsustainable O (10^ s), i.e. more than one day! Even at the 
Tevatron, where the multiplicity is quite modest, the fact that noise may cause the number 
of active calorimeter cells to be far larger than the number of particles has led to the use 
of a complex (and physically questionable) preclustering procedure prior to running the kt 
jet hnder, so as to reduce the effective value of N to something that is manageable [20]. 
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The slowness of the kt jet-finder has been one of the motivating factors behind proposals 
for alternative jet-hnders [9,10]. Here we will show that the kt jet-finder can in fact be 
formulated in an algorithmically fast (iVlniV) manner. A C++ implementation of this (and 
a related algorithm^ will be shown to run orders of magnitude faster than currently 
available implementations, making it feasible (and easy) to use the kt jet hnder for efficiently 
studying high-multiplicity events. 


2 The kt jet-finder 

The kt jet hnder, in the longitudinally invariant formulation suitable for hadron colliders, 
is dehned as follows. 

1. For each pair of particles i, j work out the kt distance dtj = mm{ktt, ktj)R‘fj with 

= iVi ~ VjY + ~ where ku, rji and 0* are the transverse momentum, 

rapidity and azimuth of particle i] for each parton i also work out the beam distance 
diB = kl- 

2. Find the minimum dmin of all the djj, diB- If c^min is a dtj merge particles i and j into 
a single particle, summing their four-momenta (alternative recombination schemes 
are possible); if it is a diB then declare particle i to be a hnal jet and remove it from 
the list. 

3. Repeat from step 1 until no particles are left. 

There exist extensions of this basic procedure, (a) where is rescaled relative to diB by a 
user-chosen factor 1/i?^ ~ 1 or (b) where clustering is stopped when all dij,diB are above 
a jet resolution threshold dcut- We here consider only the simplest version, as given above, 
but the arguments below are identical for the general case. 

Now we reconsider the above procedure, making explicit the computational overheads 
of the various steps as implemented in standard jet hnding codes [22]. 

1. Given the initial set of particles, construct a table of all the dtj, diB- 

[O (A^^) operations, done once] 

2. Scan the table to hud the minimal value dmin of the dij, diB- 

[O operations, done N times] 

3. Merge or remove the particles corresponding to dmin as appropriate. 

[(9(1) operations, done N times] 

^‘Jet-algorithm’ is often used in the literature to refer to the choice of the rules for finding a jet; here 
instead ‘algorithm’ refers to the translation of a given set of jet-finding rules into explicit steps on a 
computer. 
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4. 


Update the table of dij, diB to take into account the merging or removal, and if any 
particles are left go to step 2. 

[O (iV) operations, done N times] 


Step 2 dominates, requiring O (iV^ x N = N^) operations.^ 


3 The FastJet Algorithm 

To obtain a better algorithm we isolate the geometrical aspects of the problem, with the 
help of the following observation. 

Lemma: If i, j form the smallest dij, and ku < ktj, then Rij < for all i ^ j, i.e. j 
is the geometrical nearest neighbour of particle i. 

Proof: Suppose the Lemma is wrong and that there exists a particle i such that 
Rii < Rij'- then die = since min(A;|, k^^) < we have that da < dij, 

in contradiction with the statement that i and j have the smallest dij. 

This means that if we can identify each particle’s geometrical nearest neighbour (in 
terms of the geometrical Rij distance), then we need not construct a size-iV^ table of 
dij = mm{k^^, k^j)R^j, but only the size-iV array, dig^, where Qi is i’s geometrical nearest 
neighbour^. We can therefore write the following algorithm: 

1. For each particle i establish its nearest neighbour Qi and construct the arrays of the 
dig,^ and diB • 

2. Find the minimal value dmin of the dig^, diB- 

3. Merge or remove the particles corresponding to dmin as appropriate. 

4. Identify which particles’ nearest neighbours have changed and update the arrays of 
dig. and diB- If any particles are left go to step 2. 

This already reduces the problem to one of complexity for each particle we can hnd 
its nearest neighbour by scanning through all O {N) other particles [O {N"^) operations]; 
calculating the dj^., diB requires O (iV) operations; scanning through the dj^., diB to hnd the 
minimal value dmin takes O (iV) operations [to be repeated N times]; and after a merging 
or removal, updating the nearest neighbour information will require O {N) operations [to 
be repeated N times] 

^One notes also the storage requirement in step 1 of + O (N) bytes (double precision), which is 
manageable for N = 1000 but becomes an issue in heavy-ion environments with up to 50000 particles. 
At the (substantial) expense of recalculating the O (iV^) dij at each iteration, the storage issue can be 
eliminated. 

•^We shall drop ‘geometrical’ in the following, speaking simply of a ‘nearest neighbour’. 

■^This last point is not strictly speaking trivial: when particle i is removed or merged we have to update 
the nearest neighbour information for all particles that previously had i as their nearest neighbour — 
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We note, though, that three steps of this algorithm — initial nearest neighbour iden- 
tihcation, hnding at each iteration, and updating the nearest neighbour information 
at each iteration — bear close resemblance to problems studied in the computer science 
literature and for which efficient solutions are known: 

• Given an ensemble of vertices in a plane (specihed by the T]i and (pi of the particles), to 
find the nearest neighbour of each vertex one can use a structure known as a Voronoi 
diagram [23] or its dual, a Delaunay triangulation. The Voronoi diagram divides the 
plane into cells (one per vertex), such that every point in the cell surrounding a vertex 
i has i as its nearest vertex. The structure is useful for nearest-neighbour location 
because the vertex Qi nearest to vertex to i is always in one of the (few, i.e. O (1)) cells 
that share an edge with the cell of vertex i. An example is shown in hgure^ Voronoi 
diagrams for N points can be constructed with O (NlnN) operations (see e.g. [24]), 
and the nearest neighbour identihcation for all N points can then be performed with 
a further O {N) operations. 

• Dynamic insertion and removal of a point in the Voronoi diagram, and correspond¬ 
ing updating of all nearest neighbour information, can be performed with O (In V) 
operations [25] (to be repeated N times). 

• The array of djg. changes only by O {1) entries per iteration. Therefore one can 
represent it with a binary tree structure, whose construction requires 0{N\rLN) 
operations and in which finding the minimal value, insertion and removal are all 
guaranteed to require at most O (In N). The binary tree is constructed once at start¬ 
up, and there are then O {N) updates and searches for the minimum, leading to a 
total of O {N In N) operations. 

Therefore both the geometrical and minimum-hnding aspects of the kt jet-finder can be 
related to known problems whose solutions require O {N In N) operations. 


4 Timings 

The program Fast Jet® codes concrete implementations of the and NlnN algorithms 
described above. It has been written in C-I--I- and for the V In V case makes use of a number 
of pre-existing components. Construction of a Voronoi diagram is a sufficiently common 
task (useful in areas of science ranging from astronomy to zoology) that several codes 
are publicly available. Of these, the only one that we found that also straightforwardly 
allows the addition and removal of points from a pre-constructed Voronoi diagram, was the 

fortunately one can show that on average, the number of particles that had * as a nearest neighbour is 
of D (1). One also needs to establish if any particles acquire the newly created particle £ as their nearest 
neighbour — this can be done in O (N) time by comparing each particle’s current nearest neighbour 
distance with its distance to i. 

^Available from http;//www. Ipthe . jussieu. fr/~salam/fast jet 


5 





Figure 1: The Voronoi diagram for ten random points. The Delaunay triangulation (red) 
connecting the ten points is also shown. In this example the points 1, 4, 2, 8 and 3 are the 
‘Voronoi’ neighbours of 7, and 3 is its nearest neighbour. 



Figure 2: The running times (on a 3 GHz Pentium 4 processor with 1 GB of memory, 512 kB of 
cache, and version 3.4 of the GNU g++ compiler) of the KtJet [22] and FastJet implementations 
of the /cfclustering jet-finder versus the number of initial particles. Different values of N have 
been obtained by taking a LHC dijet event with pt — 60 GeV and adding on variable numbers of 
minimum bias events. Both kinds of events have been simulated with Pythia 6.3 [28]. 


Computational Geometry Algorithms Library (CGAL) [26], in particular its triangulation 
components [27].® For the binary tree we made use of a (red-black) balanced tree.^ 

®One issue relates to the fact that we need nearest-neighbour location on a cylinder {rj-cj) space) whereas 
CGAL works on the Euclidean plane. This problem can solved by making mirror copies of a small 
(^ 1/\/N) fraction of the vertices near the 0 — 27r border. 

^Balanced trees are the basis of the map and set classes in the C-I--I- Standard Template Library. 
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Figure 121 shows the running times for the two algorithms in Fast Jet as well as for Kt Jet, 
a standard implementation [22] of the algorithm. Our algorithm” actually departs 
slightly from exact behaviour owing to certain further optimisations carried out.^ The 
scaling with N of the Voronoi-based algorithm has been verihed to go as iVlnA^, as ex¬ 
pected. It is the fastest algorithm only for N > 10^, owing to a large coefficient in front 
of its A^lnA^ behaviour, mostly associated with the computational geometry tasks. This 
situation could conceivably be improved in the future by optimisations of the CGAL pack¬ 
age or by replacing it with a dedicated implementation of the construction and updating 
of the Voronoi diagram. 

The better of the and N In N algorithms (which can be selected based on the value 
of N) runs at least an order of magnitude faster than the algorithm for all values of N 
shown, vastly more at large N. 

Figure El compares the running time of our combined N‘^-N\nN Fast Jet implementa¬ 
tion of the kf jet-hnder with other jet-hnders whose code is publicly available. One sees 
that it runs at least an order of magnitude faster than all others (except the almost IR 
unsafe JetClu). 


5 Perspectives 

Since the Fast Jet algorithm is functionally equivalent to the standard algorithms used 
for the kt jet hnder till now, the results of the clustering are of course identical to those 
of other implementations. Howewer, its enhanced speed opens up new ways of using kt 
clustering in the analysis of high-multiplicity events. 

Historically, one apparent drawback of fct-type jets with respect to cone-type jets in 
hadron-hadron collisions was considered to be the larger fluctuations of the areas of the 
jets dehned by the clustering procedure. Such fluctuations would seem to make it more 
difficult to subtract, from the hard event, the energy coming from the non-perturbative 
underlying event and from any additional minimum bias interactions taking place in the 
same bunch crossing (pileup). 

However, the fluctuations become irrelevant if one can easily estimate the area of each 
individual jet. This can be done on an event-by-event basis, as follows: because of the 
infrared safety of the kt jet-hnder algorithm, one can add a large number of extremely soft 
particles (“ghosts”) to the event without modifying the properties of the hard jets. After 
clustering, each jet will contain a subset of the ghosts, and if the ghosts had a uniform 
density in rj and 0, then the number of ghosts in a given jet will be a measure of its area. In 
practice we have found that the use of O (10^) ghost particles is necessary to obtain reliable 
area estimations. This dehnition for the area of a kt jet can of course be implemented with 
any coding of the jet-hnder. It is however impractical, indeed nearly impossible, to deal 

®The coefficient of iV^ can be reduced by tiling the plane into rectangles of edge length > 1. Then 
for each vertex i one can limit the nearest neighbour search to its own tile and adjacent tiles — vertices 
further away will have Rij > 1 and so dij > dis ■ 
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Figure 3: The running times of various jet-finders versus the number of initial particles. JetClu 
is a widely-used cone-type jet-finder, however it is ‘almost infrared unsafe’, i.e. perturbative 
predictions have large logarithmic dependence on small parameters (e.g. seed threshold) [29,30]. 
MidPoint [29] is an infrared safe cone-type jet finder. For both we use code and parameters from 
CDF [31]. The Optimal Jet Finder [9] (OJF) has been run with ^cut = 0.15 and a maximum of 
8 jets, so as to produce a final state similar to that returned by the kt and cone jet-finders and 
to limit its run time. 


with the required large number of ghost particles without a fast code. 

Preliminary studies have shown that with simple assumptions about the uniformity of 
the underlying event and pileup, one can readily determine its size and subtract it from the 
hard jets, leading to good determinations of kinematical quantities (e.g. invariant masses) 
in high-luminosity pp collisions, or of single inclusive jet distributions in Pb-Pb collisions 
at the LHC. Full results will be shown in [32]. 

Two more observations are worth making before closing this section. They will both 
be discussed in more detail in [32]. 

The first is that it can also be interesting to examine alternative definitions of jet 
areas. One option is to make use of the areas of the Voronoi cells of all the real particles 
belonging to a given jet. This definition avoids the need to cluster thousands of ghost 
particles together with the real ones. It instead rests on the geometrical properties of the 
event, and on the computational geometry component of the Fast Jet implementation. 

The second observation is that there exist clustering-type jet-finders other than the kt 
jet-finder that share a large fraction of its features (including, of course, infrared safety), 
and the possibility of a fast implementation. One such example is the “Cambridge” jet- 



finder. It was originally formulated in the context of collisions in [33] and an inclusive 
version, adapted to hadron collisions, was given in [34], We shall call this inclusive version 
the Cambridge/Aachen algorithm. It is dehned in the same way as the kt jet-hnder at the 
beginning of Section El but replacing the particle-particle distance by dij = and 

the particle-beam distance by diB = 1- We shall show in [32] that the Cambridge/Aachen 
jet-hnder has smaller average areas than the kt jet-hnder, making it perhaps even better 
suited for jet studies in high-multiplicity environments. 


6 Conclusions 

To conclude, we have identihed an unexpected relation between clustering type jet-hnders 
and a widely studied problem in computational geometry. The resulting reduction of the 
complexity of the kt jet-hnding problem, from to A^lnA^, opens up the previously 
inconceivable option of using the kt jet-hnder for the large values of N that arise when 
considering all cells of a hnely segmented calorimeter and for heavy-ion events. For mod¬ 
erate N, the one or two orders of magnitude in speed that we gain with a related 
approach pave the way to much wider use of the kt jet hnder for normal hadron-collider jet 
studies, especially at the LHC. More generally, the speed gains discussed in this paper also 
suggest novel ways of using the kt jet hnder, which are the subject of ongoing investigation. 
One example, given in section El is the determination of jet areas, knowledge of which is 
crucial for optimal subtraction of pileup contamination in high luminosity environments. 
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